The shearing instabiUty of a dilute granular mixture 

J. Javier Brey and M.J. Ruiz-Montero 

Fisica Teorica, Universidad de Sevilla, 
Apartado de Correos 1065, E-4IO8O, Sevilla, Spain 
(Dated: February 20, 2013) 

Abstract 

The shearing instability of a dilute granular mixture composed of smooth inelastic hard spheres 
or disks is investigated. By using the Navier-Stokes hydrodynamic equations, it is shown that 
the scaled transversal velocity mode exhibits a divergent behaviour, similarly to what happens in 
one-component systems. The theoretical prediction for the critical size is compared with direct 
Monte Carlo simulations of the Boltzmann equations describing the system, and a good agreement 
is found. The total energy fluctuations in the vicinity of the transition are shown to scale with 
the second moment of the distribution. The scaling distribution function is the same as found in 
other equilibrium and non-equilibrium phase transitions, suggesting the existence of some kind of 
universality. 

PACS numbers: 45.70.-n,05.20.Dd, 05.60.-k,51.10.+y 
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I. INTRODUCTION 



The simplest statistical mechanics model for granular gases is a system of smooth inelastic 
hard spheres or disks, with constant coefficient of normal restitution [l, 2]. One characteristic 
feature of these systems, as compared with their elastic limit, is their instability against 
small wave-vector spatial fluctuations when they evolve freely [Jlsl. In the case of one- 
component granular gases, this spontaneous symmetry breaking leads to the formation of 
velocity vortices and density clusters, being often referred to as the shearing or clustering 
instability of the homogeneous cooling state (HCS). It is accurately predicted by a linear 
stability analysis of the hydrodynamic Navier-Stokes equations of granular gases, which 
also shows that it is driven by the transversal shear mode Moreover, fluctuating 

hydrodynamics is able to predict not only the initial set-up of the spatial correlations {7,8], 
but also the behaviour of the system near the critical point of the instability js]. This 
includes the critical exponents governing the behaviour of both macroscopic properties and 
fluctuations. 

In the last years, the hydrodynamic theory of granular gases has been extended to binary 
mixtures. Navier-Stokes equations for the hydrodynamic fields of the mixture, with explicit 
expressions for the involved transport coefficients, have been derived 1^, 11 1. The hydro- 



dynamic equations for a mixture are much more involved than those for a one-component 
system and, therefore, so is the hydrodynamic linear stability analysis of the HCS 12|. But 
there is no reason to expect that the physical mechanisms leading to the shearing instability 
in simple granular gases does not hold for mixtures. If that is the case, the behavior of the 
transversal component of the velocity field as the size of the system increases is the origin 
of the instability. And it happens that the evolution equation for this hydrodynamic mode 
is decoupled from the equations for the rest. This feature greatly simplifies the analysis of 
the initial set up of the instability. 

The aim of this paper is to investigate the behavior of the transversal velocity mode in a 
binary granular mixture identifying, in particular, the existence of the shearing instability 
and determining the critical point predicted by the hydrodynamic theory. Also the behaviour 
of some average properties of the granular mixture in the vicinity of the instability will be 
studied and compared with those of a single component granular gas. In addition, the 
critical behaviour of some quantities that are peculiar of granular mixtures, as the non- 
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equipartition of kinetic theory will be investigated. It is worth to emphasize that both 
theory and simulations presented here are restricted to the linear regime in which deviations 
of the fields from their values in the HCS are small. 

The remaining of this paper is organized as follows. In Sec. HIl the relevant properties 
of the HCS of a granular mixture in the context of kinetic theory are summarized. This 
includes the criterion determining the partial temperatures of both components of the mix- 
ture. Section IIIII consists of a short review of the linear stability analysis of the transversal 
velocity field of the HCS to Navier-Stokes order [l2|. The associated eigenvalue is identified, 
and it is shown that it has a qualitative change of behaviour when the size of the system is 
larger than a critical value, which depends on the parameters defining the system. In Sec. 
IIVI the dynamics of the inelastic hard spheres or disks is reformulated by means of a change 
of variable, so that the HCS is mapped onto a steady state. It is a straightforward exten- 
sion of a method previously developed for one-component granular systems. This steady 
representation is used in Sec. |V] to perform direct Monte Carlo simulations (DSMC) of the 
system, whose results are compared with the theoretical predictions. The behaviour of the 
total energy fluctuations is also analyzed using the simulation results. It is seen that the 
relative dispersion of the energy fluctuations exhibits a power-law divergent behaviour near 
the instability. The paper concludes with a short discussion of the results and an analysis 
of the shape of the probability density distribution for the total energy fluctuations. 



II. THE HOMOGENEOUS COOLING STATE OF A DILUTE GRANULAR MIX- 
TURE 

A fluidized binary mixture of smooth inelastic hard spheres {d = 3) or disks {d = 2) 
is considered. The mass and diameter of particles of species i {i = 1,2) are mj and (Tj, 
respectively. The inelasticity of collisions is assumed to be described by constant, velocity 
independent, coefficients of normal restitution. There are three of them: an, 022, and 
ai2 = 0:21, where aij refers to the collision of a particle of species i and a particle of species 
j. These coefficients are defined in the interval < aj-,- < 1, the value unity being the limit 
of elastic collisions. 

The macroscopic fields number densities ni{r,t), flow velocity u{r,t), and granular tem- 
perature T(r, t) are deflned in the usual way as local velocity moments of the distribution 
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function of the system, although setting the Boltzmann constant equal to unity. More pre- 
cisely, they can be expressed in terms of the one-particle distribution functions of the two 
species fi{r,v,t) as 

niir,t) = [ dvfi{r,v,t), (1) 

(2) 



n 



p{r,t)u{r,t) = ^ / dvmivfi{r,v,t), 

i=l,2 
i=l,2 ^ 



(3) 



where p = mini + 7712^2 is the total mass density, n = ni -|- n2 is the total number density, 
and V = V — u is the peculiar velocity. 

In this paper, attention will be restricted to a low density system. Then, the time 
evolution of the one-particle distribution functions is given by a pair of coupled nonlinear 
Boltzmann equations 

{dt + v- V) /,(r, v,t)=J2 Ah. /,] , (4) 

i=i,2 

i = 1,2, and Jij denoting the Boltzmann collision operator describing the scattering of pairs 
of particles From Eqs. (jlj), balance equations for the macroscopic fields are derived by 
multiplying by 1, v, and v'^, respectively, and subsequent integration over the velocity. They 
have the form 

dtrii + V-{nu + ji) = 0, (5) 

dtU + u-Vu + p~^V - P = 0, (6) 

dtT + u-VT--V -J" ji + ^iV -q + P :Vu)+TC = 0. (7) 
n ^-^ nd 

i 

In the above expressions, ji is the number of particles flux for species i relative to the 
local flow, P is the pressure tensor, q is the total heat flux, and Q is the cooling rate giving 
account of the loss of energy in collisions. These quantities are defined as functionals of the 
one-particle distribution functions. 

The balance equations ([5])- ([7]) admit time- dependent homogeneous solutions characterized 
by uniform densities nj /j, a vanishing velocity field Uh = 0, and an homogeneous granular 
temperature evolving in time accordingly with 



dtTh{t) = -a{t)Th{t), 
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(8) 



where (h is the coohng rate of the homogeneous state. Of course, this equation is only 
meaningful if the cooling rate is expressed in terms of Uh and Th, then becoming a closed 
equation for T^. This is accomplished when the distribution functions of both species depend 
on time only through the granular temperature Th{t). The distribution functions fi h{v,t) 
having this property are said to be "normal " and define the homogeneous cooling state 
(HCS) of the mixture [l3j, which is the state considered in this paper. 

Partial temperatures of the species in the HCS, Ti^h{t), are defined through 



ni,hTi,h(t) 



Therefore, it is 



ni,hTi,h{t) = nhTh{t). 



(9) 



(10) 



i=l,2 



Evolution equations for the partial temperatures are directly derived from the Boltzmann 
equations, particularized for the HCS, 

dtTi^h{t) = ~Ci,h{t)Ti,h{t), (11) 
with the partial cooling rates Ci,h{t) given by 

Ci,h{'t) = r} / dvmiv'^Jij[v\fi^hJjA^ 

ni,h-Li,h{t)a ^—^ J 

3 

Consistency of Eqs. ([8]), f fTOj) . and f lTTj) requires that 

nhTh(t)(h{t) = ni^hTi,h{t)C,i,h{t)- 

1=1,2 

Moreover, as a consequence of the distribution functions of the HCS being normal it is [1 



(12) 



(13) 



CiAt) = C2,h{t) = Chit). 



(14) 



The explicit evaluation of the cooling rates requires us to solve the coupled pair of Boltz- 
mann equations for the distribution functions of the HCS. Nevertheless, a quite accurate 
approximation, at least for not very strong inelasticities, is obtained by using Gaussian dis- 
tributions for fi h{v,t) with the second moments corresponding to the partial temperature 
Ti^hit). In this way, it is obtained |lOl. Il3l-ll5|. 



i=i,2 



9i9j 



1/2 



[l+a 



1 - ^ij. 



l + a. 



+ 



9. 



(^12 



d-1 



(15) 



where Xi = rii/n is the number concentration of species i, aij = ((Xi + crj)/2, 

vo{T)^[—] , 



with 



is a thermal velocity, 



is a characteristic length, 



and 



fx ) 

mxvn^ 
mi + 1712 ' 

-l\l/2 



nrti 



rrii + nij 



(16) 

(17) 
(18) 
(19) 
(20) 



Now, the expressions of and (2,h can be introduced into Eq. (fT4|) . The solution of 
the resulting equation provides 7 = T-\h{t) jTo hit) as a function of ni^h/n2^h and the other 



parameters defining the model [lO|, ll3|, llSl, ll6j. The accuracy of the results derived in this 
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19| . as well as 



way has been checked by comparing with molecular dynamics results 
with data obtained from particle simulations of the Boltzmann equations using the DSMC 
method. 



III. LINEAR ANALYSIS OF THE TRANSVERSAL SHEAR MODE OF A GRAN- 
ULAR MIXTURE 

In the case of a one-component granular system, it is well known that the HCS becomes 
unstable when a linear size of the system is larger than the critical value, which depends on 
the parameters defining the state 31,^]. For dilute granular gases, a linear stability analysis 
of the hydrodynamic Navier-Stokes equations shows that the origin of the instability lies in 
the behaviour of the transversal velocity mode. For sizes of the system larger than the critical 
one, the transversal velocity decays slower than the square root of the temperature {2^. As 
a consequence, nonlinear coupling of the scaled modes becomes relevant and a clustering 
instability develops in the system. It is worth to stress that the transversal component of 
the velocity field is not itself linearly unstable, but it enslaves the other hydrodynamic modes 
through some nonlinear coupling. In particular, this leads to an increase of the temperature 



in the regions of larger vorticity. Then, a pressure gradient shows up and produces a density 
fluctuation leading to the formation of clusters. A detailed account of the theory and the 
comparison with simulation results is given in {g]. 

For a binary granular mixture, the complexity of the Navier-Stokes hydrodynamic equa- 
tions 
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llj makes it more complex the full linear stability analysis of the HCS. Never- 



theless, in the mixture the scaled transversal velocity mode is decoupled from the other 
hydrodynamic modes, as it happens for one-component granular gases [l2^. Consequently, 
the analysis of the transversal mode is quite simple. Here the change in its time behaviour 
will be investigated. This change indicates the need of considering nonlinear couplings be- 
tween modes and can be considered as the precursor of a clustering instability, by analogy 
with the one-component case. 

Consider the evolution equation for the velocity, derived from the momentum conserva- 
tion, Eq. (Q. To Navier-Stokes order, the pressure tensor P has been computed by using the 
Chapman-Enskog method to solve the pair of coupled Boltzmann equations of the mixture. 



It reads 
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2l| 



P = p\ — rj 



(Vw) + (Vn) 



d 



y-u)\ 



(21) 



In this expression, p = nT is the local pressure, I is the unit tensor of dimension d, and r] is 
the coefficient of shear viscosity of the mixture, that can be written as 

pX 



V 



VoiT) 



V 



(22) 



where A = (na 



d-l\-l 
12 J 



and Vq{T) has been defined in Eq. (fT6l) . Moreover, rj* is a dimensionless 



function of the partial temperatures and the concentrations of the species. Its explicit form 
is given in the Appendix. 

The Navier-Stokes equation for the velocity resulting after substituting Eq. ( |2T1) into Eq. 
(|6]) will be now linearized around the HCS. Firstly, the hydrodynamic fields are written in 
terms of their deviations from the HCS values, 



u{r, t) = 5u{r, t), 
T{r,t) = n{t) + 6T{r,t). 



(23) 

(24) 
(25) 



7 



Keeping only up to first order in the deviations of tlie fields, the Navier-Stokes equation for 
the velocity becomes 



Ph 



dSu 



0, 



(26) 



with Ph and rjh being the mass density and the shear viscosity of the reference HCS, respec- 
tively. 

At this point, it is convenient to introduce dimensionless time and space coordinates 
such that the time dependence of the reference state be eliminated in Eq. (|6]). Then, new 
variables are defined by 

(27) 



In the time scale r, the evolution equation for the temperature of the HCS becomes 

drTh{r) = -Cnir), 
with the time-independent reduced cooling rate (* given by 

Therefore, the temperature of the HCS decays exponentially on the r scale. 
Also, reduced hydrodynamic fields are introduced, 

Su 



Sn 

rih Vo{Th 



Use of the definitions in Eqs. f l27|) and ( 130|) into Eq. (!26|) yields 

2 



with 



/9 = 



dy d-2 d fduj 
di) d di [~dr 



2{xiTni + X2m2) 
Now, the Fourier representation defined as 



h = J dle-^'^-'fil), 
for an arbitrary function f{l), will be employed. The transformed of Eq. fl3ip is 



(28) 



(29) 



(30) 



(31) 



(32) 



(33) 



(34) 
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Consider the transversal component uJk^± of the velocity field, i.e. the vector component of 
ujk perpendicular to k. Its evolution equation is trivially obtained from Eq. f lM|) and has 
the form 

%i+(w*^-f)a.,x = 0.. (35) 
This is a closed equation for uJk,±, similar to that found in one-component granular gases 



20j |. An equivalent equation has been derived in [12|, where the general issue of the linear 
stability analysis of the Navier-Stokes equations of a dilute granular gas is addressed. The 
solution of Eq. f l35|l reads 

= e-^^^c:;fc,±(0) (36) 

where the decay rate s± is 



^ (5r^*e - ^ . (37) 
This leads to the identification of a critical value of the wavenumber vector given by 

/ A* \ 1/2 

kc=Uh:] ■ (38) 



A linear excitation of the scaled transversal velocity with k < kc grows in time. Therefore, 
vortices of the scaled velocity field are expected to develop in time when excitations of this 
kind are present in the system. This does not mean that the actual velocity field u is 
linearly unstable. In fact, it is easily seen that the perturbation 6u decays exponentially in 
time because of Eq. ( I28l) . The result above indicates that the linear analysis will eventually 
fail and nonlinear effects associated to coupling of hydrodynamic modes will have to be 
taken into account. In simple granular gases, this coupling leads to the development of the 
clustering instability. For this reason, the formation of vortices in the scaled velocity field 
is sometimes referred to as the shearing instability of the HCS. 

IV. MAPPING THE HCS ONTO A STEADY STATE 

In order to verify the validity of the ideas developed above and, in particular, to check 
whether the shear instability also exists in granular mixtures and if it is accurately predicted 



by the hydrodynamic Navier-Stokes equations, the DSMC method |22l-l24| has been used 
to generate numerical solutions of the coupled pair of Boltzmann equations. Actually, the 
method is an A^-particle algorithm designed to mimic the dynamics of a low density gas and, 
therefore, it also provides equilibrium and non-equilibrium fiuctuations and correlations. 



One of the technical advantages of the DSMC method is that it permits to incorporate 
in the simulations the symmetries of the particular situation of interest. This allows a sig- 
nificant increase in the statistical accuracy of the measured properties. Here the aim is to 
investigate the development of inhomogeneities in the vicinity of the critical size associated 
with the shearing instability. Therefore, the simulation must allow the formation of sponta- 
neous fluctuations of a given wavelength. For this reason, it is enough to consider a system 
in which gradients can occur in only one direction, arbitrarily taken as the x axis. The 
components of the position of the particles perpendicular to that axis are not relevant from 
the point of view of the simulation. In other words, the simulation is restricted to systems 
which are homogeneous in the planes perpendicular to the x axis. The system size in the x 
direction is L, and periodic boundary conditions are used in that direction. 

A simulation of the cooling mixture in the actual phase space variables is difficult since the 
rapid cooling of the system leads to rather small energies, and large uncertainties very soon. 



26| for one- component granular gases 



To deal with this, the procedure introduced in [25 

and extended to mixtures in 1^, will be employed. The idea is to exploit the existence of 
an exact mapping of the HCS onto a steady state. Although the method can be formulated 
in the time scale r defined in Eq. fj271) . this would have the technical complication that the 
exact cooling rate is not known a priori. Consequently, it is convenient to introduce a new 
time scale s by 

ds=^dr= ^-^^ dt, (39) 

where w is an arbitrary dimensionless frequency. Now, the positions and velocities of the 
particles are represented in the I and s scales. The particle dynamics in these variables 
consists of an accelerating streaming between collisions, 

I - («) 

while the effect of the collision of two particles is the same as in the original time scale, 
given its instantaneous character. The dynamics defined by Eqs. ( HOl) and ( HTj) is seen to be 
equivalent to a change in the original time scale, 

t^7s = ln-, (42) 

^0 
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where to is another arbitrary constant. The acceleration term in the dynamics ( BTj) is able 



to balance the energy lost in collisions, thus enab 
temperatures in the new dynamics are given by [19 1 



ing a steady state. The steady partial 



(43) 



From Eqs. ffT^ and ffUl) . it follows that also 



o \ 2 



t: = J ' (44) 

= CiT) /T^/"^ . The above mapping does not affect the hydrodynamic shear instability. In 
the time scale s, Eq. (|35|) becomes 

+ t: ujk± = 0. (45) 



ds \ C* 

As expected from dimensional analysis, the arbitrary constant w plays no role in the stability 
criterion. 

Because of the symmetry of our simulations as described above, it is clear that the 
minimum wavevector allowed is given by kmin = ^irXh/L. Therefore, the stability condition 
given by Eq. ( l38l) is equivalent to L < Lc with the critical length Lc given 

by 

L^ = 2nxJ^y^\ (46) 



It is worth to mention the existence of a related instability for the total momentum of the 
system in the scaled variables Nevertheless, it is not physically relevant and can be 

eliminated by taking a vanishing initial total momentum. 



V. SIMULATION RESULTS 



In the simulations to be reported in the following, a binary mixture of = A^i + A''2 
inelastic hard spheres {d = 3) has been used. In order to reduce the number of parameters 
characterizing the system and to allow for a systematic study of those being varied, the 
number of particles of both species, and therefore the concentrations, have always been the 
same (A"! = A'2), as well as the diameters of the particles, i.e. cxi = (J2. Moreover, the 
coefficient of normal restitution for collisions between particles of different species has been 
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taken as the average of the coefficients for equal species colhsions, i.e. ai2 = (an + a;22)/2. 
On the other hand, the mass ratio A = 1712 /mi and the two coefficients of normal restitution 
ttii and a22i as well as the size L of the system have been varied in the simulations. 

The behaviour of several properties of the HCS as the size L of the system approaches the 
critical value has been studied. The number of particles per unit of length in the x direction 
has been kept fixed in the simulations, = N/L^ = 2000£^^, where ih = Xh/V^t^ = 
(v^TTCT^n/j) ^ is the mean free path. It is worth to stress that the number of particles used 
in the DSMC method does not affect the validity of the low density limit, that is inherent 



to the method itself 



'221, 



The simulations were performed using the steady representation of the HCS discussed in 
the previous section and, unless explicitly otherwise established, the values of the properties 
reported in the following have been averaged in time once the system had reached the 
steady state, as well as over a number of different trajectories (typically 50). The value of 
the arbitrary constant w was chosen in all cases as w = Cg/2, with (q being the value of 
C obtained in the Gaussian approximation, i.e. those given by Eq. ( IT5l) . If the Gaussian 
approximation were exact, the measured steady value of the total temperature would have 
been one. 

The ffist point addressed in the simulations was to check that the scaled transversal 
velocity field was really the ffist hydrodynamic mode becoming unstable as L increases. 
The steady state reached by the system for different sizes was investigated, starting from 
a system with L small enough as to guarantee that the HCS was stable, and increasing L 
from there on. The different hydrodynamic fields were monitored at different times, and in 
all cases it was found that they were the y and z components of the scaled velocity field 
the ffist hydrodynamic modes exhibiting large fluctuations. As long as L is not large, these 
fluctuations eventually decay, but when the system size was increased enough, a non-decaying 
scaled transversal velocity field emerged. This is illustrated in Fig. [H where the density and 
one of the components of the transversal velocity field are shown at three different times 
for a system with an = 0.8, ^22 = 0.98, m2/mi = 4, and L = 23.7ih- Both the total 
hydrodynamic fields as well as those associated with each of the species are displayed. The 
velocity field for each of the components are defined by equations similar to Eq. (|2]) j^^], 
and they have been scaled with the square root of the temperature of the system. It is seen 
that the system exhibits an spontaneous perturbation of the transversal velocity field that 
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FIG. 1: (Color online) Snapshots of the dimensionless scaled density (left) and of the one of 
the components of the perpendicular velocity field ojy for a system with an = 0.8, 022 = 0.98, 
m2/rni = 4, and L = 23.7ih- The (blue) circles correspond to the whole fluid properties, while the 
empty (black) square and (red) triangles are for species 1 and 2, respectively. From top to bottom 
the times are s = 6.68 x 10'^, 1.33 x 10^, and 1.99 x 10^, in the dimensionless scale defined by Eq. 
(j39|) . with zu chosen as discussed in the main text. 



does not decay in time, and corresponds to the first possible harmonic. On the other hand, 
the density remains homogeneous. Note that the local average velocities of the species is 
the same as that of the whole fluid. Then, it was concluded that the scaled shear mode is 
the field for which the linear approximation first breaks down. 

To measure the critical size Lc, the following procedure was used. First, the average 
value of the total energy in the steady state E was measured as a function of the size of the 
system. Then, it was assumed, to be checked in the simulation results, that the behaviour 
of the average steady energy near but below the shearing instability obeys a law of the form 



< E > - < E >h 
oe = ^ oc 

< E >h 



Lc-L 

Lr 



(47) 
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FIG. 2: (Color online) Relative deviations 6e = {< E > — < E >h)/ < E >h of the average total 
energy of the system from its asymptotic values in the HCS, as a function of the system size L, in 
the vicinity of the shear instability. The parameters of the system are an = 0.92, a22 = 0.98, and 
1112/ mi = 8. The symbols are from the simulations and the straight lines are fits in the "critical 
region". The (blue) circles correspond to the whole fluid, while the empty (black) squares and 
(red) triangles are for species 1 and 2, respectively. 



where < f >h denotes the constant asymptotic average value of the property / in the HCS, 
far away from the shear instabiUty, also obtained from the simulations. The above behaviour 
is suggested by the results obtained for a one-component dilute granular gas near its shear 
instability |28|,|29|. 



In Fig. [21 5^ is plotted as a function of the system size. The observed linear behaviour 
is consistent with Eq. f l47p . and from the parameters of the linear fits, simulation values of 
the critical size Lc are directly obtained. The fits of the three lines lead to the same value, 
namely — 30.13ih- A similar behaviour was obtained for all the values of the restitution 
coefficients and the mass ratio investigated. It follows that the increase of the total average 
energy of the system and also that of each of the components, as the system approaches the 
instability, is characterized by Eq. (H71) . 

The comparison between the measured critical sizes and the theoretical prediction given 
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FIG. 3: (Color online) Dimensionless critical size Lc/d-h for the shear instability of the HCS as a 
function of the mass ratio 7712/7711. The symbols are from the DSMC method, while the solid lines 
are the theoretical predictions given by Eq. (j46|) . Three different sets of values of the coefficients 
of normal restitution have been considered, as indicated in the insets. 

by Eq. ( H6i) is presented in Fig. [3] as a function of the mass ratio 777,2/^1, for three different 
sets of values of the coefficients of normal restitution. The agreement is quite good over the 
two decades considered. The non-monotonic dependence of the critical length on the mass 
ratio for given coefficients of normal restitution must be noticed. This is specially relevant 
for strong inelasticities. 
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FIG. 4: (Color online) Temperature ratio 721 = T2/T1 in the steady state divided by its HCS value 
721, /i) as a function of the relative "distance" to the critical size, 6l- The coefficients of normal 
restitution are an = 0.92 and 022 = 0.98. Results for different mass ratios A = m2/mi are shown, 
as indicated in the inset. 

Another quantity investigated in the simulations is the temperature ratio, 721 = T2/T1. 
Notice that as a consequence of Eqs. (IT^ and (H3|) . it is T2^h(t)/Ti^h(t) = T2g/T^g. For 
mixtures whose components have very dissimilar masses, a small but systematic deviation 
of 721 from its HCS value, 72i./t, was observed when the system approaches its critical size. 
This deviation is larger the closer the length of the system to its critical value. It is found 
that 721 > 721, h for m2 > mi, while 721 < '~i2i,h for m2 < mi. Finally, for equal masses of 
both components, no deviation from the HCS value is observed. In any case, it must be 
noticed that the deviations from the HCS values are never larger than 1%. This behaviour is 
shown in Fig. HJ in which 721 is plotted as a function of the reduced distance to the critical 
point 5l = {Lc — L)/Lc for a system with an = 0.92 and 022 = 0.98. Results for several 
values of the mass ratio, A = 1712/ mi are displayed, as indicated in the figure caption. 
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VI. FINAL COMMENTS 



The results presented in this paper show that a freely evolving dilute granular mixture 
exhibits an instability associated to the transversal shear modes that is similar to the one 
occurring in one-component granular gases. The existence of the instability, and the parame- 
ters characterizing the critical point, are accurately predicted by the linearized hydrodynamic 
equations to Navier-Stokes order. Although the detailed nonlinear mechanisms leading to 
the formation of density clusters beyond the shear instability have not been investigated 



here, it seems evident that they are the same as those for one-component systems 



a], given 



the similarity of the behavior of both, mixtures and simple systems, when approaching the 



critical size 



3Q|. 



Another relevant quantity to characterize the system near the instability is the second 
moment of the fluctuations of the total energy, defined as 

The factor has been introduced to scale out the dependence due to the number of particles 
(or size L) of the system 28|. Consider 



- (49) 

where S/j is the steady value of S far away from the instability. In one-component gases, it 
was found that near but below the shearing instability, 

(5s2 oc 5f, (50) 

with Sl given by Eq. (H7I) . Our simulation results clearly indicate that this relation also 
holds for mixtures. Actually, if the numerical results for the energy dispersion are fitted to 
it, and the fitting parameters are used to determine the critical length L^, the values are 
the same, within the statistical errors, as those obtained from the critical behaviour of the 
average energy and discussed above. 

To analyze in more detail the nature of the energy fluctuations, its distribution function 
was measured in the simulations. Again prompted by the properties of the critical region in 



one-dimensional granular gases |28|, the quantity 



E-<E> 

^~ <{E-<E >)2 >i/2 ^^-"^ 
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is considered. Far away from the instability, i.e. L <^ Lc, the probabihty distribution of 
e is Gaussian, as expected. Nevertheless, as the instability is approached, the distribution 
strongly deviates from a Gaussian, showing a clear asymmetry around the mean value. 
Moreover, and quite surprisingly, close enough to the critical length, the data for different 
values of the restitution coefficients, of the mass ratio, and of the length of the system, 
collapse onto the same curve, as illustrated in Fig. |5l This indicates that all the dependence 
of the distribution on the parameters of the system occurs through the second moment of 
the distribution. In addition, the shape of the distribution is the same as that found in one- 



component granular gases 



31 



as well as in other equilibrium and non-equilibrium systems 



32| . In these cases, an accurate expression to fit the data is 

Po{y) = K{e^-^y'\ x = h{y-s). (52) 

The values of the parameters in the above distribution follow from the normalization, zero 
mean, and unit variance conditions, with the result K = 2.14, b = 0.938, and s = 0.374. 
Therefore, it has no fitting parameters. In Fig. [5] the solid line is the plot of PQ{—e). A 
remarkable agreement with the simulation data is obtained. A peculiarity of the present 
case as compared with the other systems in which the distribution fl52|) has been used, is 
that here it is the symmetric with respect to the origin the one fitting the numerical data. 
While in granular gases, large positive fluctuations are more frequent than their symmetric. 



31 



32|. It is 



it happens the other way around in the molecular systems considered in 
possible that this difference be due to the dissipative character of granular systems. 

It is worth mentioning that in the case of one-component granular gases, the critical 
behaviour of the system near the shearing instability can be, at least qualitatively, under- 
stood in terms of nonlinear fluctuating hydrodynamics couplings [9]. It is expected that the 
analysis can be extended to binary mixtures with the same degree of accuracy. 
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0,001 



FIG. 5: (Color online) Probability density function of the scaled relative energy fluctuations e close 
to the critical size Lc- Results for several values of the coefficients of normal restitution aij, the 
mass ratio A = 1712 /mi, and the scaled length 6l = (L — Lc)/Lc are shown, as indicated in the 
inset. The symbols are from simulations, and the solid line is the distribution function given in 
Eq. ([52]) . changing y into — e. 

Appendix A: Reduced viscosity of a dilute binary mixture 



The expression for the reduced viscosity rj* of a binary mixture of inelastic hard particles 



has been obtained in 



lOj. It is given here for the sake of completeness and also because 



there is an error in the results reported in the aforementioned reference. The expression for 
the reduced viscosity reads 

V* = xi-flvl + ^2llrh, (Al) 

with 



^1 



72(2x22 - C*) - 271^2 



7i72 C*^ - 2C*(ni + T-22) + 4(riir22 - Tx2r2\\ 
The coefficients t\\ and t\2 are given by 

1/2 



2vr('^-i)/2 



Til 



CTl 



xi(2ei)~^/2(3 + 2rf - 3aii)(l + an) 



d(d+2)r(f) 1 V^i2, 

+2X2;U2l(l + a^2)QTQ2^^'' + 3)(/il2^2 " /^21^l)^r' (^1 + ^2)"'/' 



(A2) 
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ri2 = — — :::z77titX2 — Oi O2 (1 + "12) 

' '112 



rf(c/+2)r(f) V 
+ 



2d{d+l)-A , , .^-l/2 

2{d-l) ^^'^^'^ 

{d + 3)(^12^2 - 1^2161)6^^6, + 62)-'^^ 



(A3) 



3 + 2d — 3ai2 



H2l62-\6, + ^^2)^/^ - ^^^7^\^6,\6, + ^^2)-^/^ 



(A4) 



2 "'^'"^ ' 2{d-l) 
The quantities /ijj, 614, and (J12 have been defined in Eqs. f|T9l) . fl20|) . and above Eq. f|T6l) . 
respectively. Moreover, 



0"12 



(Ti + (72 



(A5) 



The expression for the reduced contribution viscosity rj^ can bo obtained from Eqs. flAip - 
l5|) just interchanging the indexes 1 and 2. 
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